Evolving momentum-projected densities in billiards 
with quantum states 
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The classical Liouville density on the constant energy surface reveals a number of interesting 
features when the initial density has no directional preference. It has been shown (Physical Review 
Letters, 93, 204102 (2004)) that the eigenvalues and eigenfunctions of the momentum-projected 
density evolution operator have a correspondence with the quantum Neumann energy eigenstates in 
billiard systems. While the classical eigenfunctions are well approximated by the quantum Neumann 
eigenfunctions, the classical eigenvalues are of the form {f{VEnVt)} where {E„} are close to the 
quantum Neumann eigenvalues and v is the speed of the classical particle. Despite the approximate 
nature of the correspondence, we demonstrate here that the exact quantum Neumann eigenstates 
can be used to expand and evolve an arbitrary classical density on the energy surface projected on 
to the configuration space. For the rectangular and stadium billiards, results are compared with the 
actual evolution of the density using classical trajectories. 

PACS numbers: 05.45.Mt, 03.65.Sq,05.45.Ac, 31.15.Gy 



I. INTRODUCTION 

The phase space density provides an alternate route for 
the study of classical dynamics that is both interesting 
and practical. Its evolution is governed by the Perron- 
Frobenius (PF) operator, £*, 



p{x, t)^C* o p{x) = S{x- f{x'))p{x')dx' (1) 



where p{x) refers to the phase space density at time t = 
0, cc = [q^p) is a point in phase space and P{x) is 
its position at time t. Equivalently, one may solve the 
Liouville equation 



(2) 



to determine the evolution of the density. A knowledge 
of the spectral decomposition of £* allows one to evaluate 
correlations, averages and other quantities of interest (T], 

For integrable Hamiltonian dynamics, the eigenvalues 
and eigenfunctions of are well known. In terms of the 
action-angle variables {1,6), the eigenfunctions are 



{2n 



,//2 



S{I-I') € 



(3) 



where /' labels the torus on which the eigenfunction has 
its support, {n\n G Z} is a measure of the spatial varia- 
tion of the eigenfunction and / is the number of degrees 
of freedom. The corresponding eigenvalues are 



(4) 



where a; (J) = dH/dl are the frequencies on the torus. 
Thus, the spectrum of an integrable system is discrete on 



a given torus /'. Note that a; generally depends on I' 
so that on the constant energy surface, the spectrum is 
continuous [l,ll[E[i|. 

The density, projected on to the configuration space, 
is the object of study in this paper. We shall be inter- 
ested in the special case where the phase space density 
is constrained to the constant energy surface. Further, 
we shall assume that the initial density has no directional 
preference and its dependence on momentum arises solely 
through its magnitude in the energy conservation require- 
ment. Such an initial density can be expressed as : 



Pe{p, q) = S{Eo - H{p, q))g{q). 



(5) 



Note that the hamiltonian nature of the flow ensures that 
the density remains confined to the constant energy sur- 
face for all times. Its projection to the configuration 
space can be expressed as 



pp(g) 



dp pE{p,q) 



(6) 



where the subscripts P and E denote 'momentum pro- 
jection' and restriction to the constant energy surface 
respectively. 

We are interested here in the time evolution of the 
momentum-projected density pp{q). Thus, we wish to 



determine the kernel Kp{q,q' 
which evolves pp{q): 



, t) of the operator C\ 



£*pOpp{q) ^ J dq' Kp{q,q',t) pp{q') 



(7) 



A straightforward way of arriving at the form of 
Kp{q, q' , t) is by projecting the initial density pE to the 
configuration space after evolving it for a time t i.e. 
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C*popp{q) = J dp o pEiq,p) 

= JdpJ dq'dp' 5{x - f{x')) 

x6{Eo-H{p',q'))g{q') (8) 

For standard hamiltonians H{p, q) = + V{q) (as- 
suming for convenience 2m = 1), the projected density 
Pp{q) = I dppE{p,q) = 7rg(q). Thus 



C'popp{q)^ J dq'Kp{q,q\t-Eo)pp{q') 



(9) 



where 



Kp{q,q',t;Eo) = ^J dpdp'S{x 



f{x'))5{Eo-H{x')) 
(10) 

and x' ~ (p', q'). Note that the kernel Kp is parameter- 
ized by Eq, the energy at which the initial phase space 
density is constrained. In general, if the starting point is 
a momentum projected density pp{q), C}p evolves pp{q) 
isotropically at an energy Eq. 

We are interested here in the time evolution of such a 
projected-density pp(q) using the eigenvalues and eigen- 
functions of £p [7|. For billiards [i, Q, the eigenfunc- 
tions of £p are well approximated by the quantum Neu- 
mann eigenfunctions, while its eigenvalues are of the form 
{fiVEnVt)} where {En} are close to the quantum Neu- 
mann eigenvalues, v being the speed of the classical par- 
ticle. Despite the approximate nature of the correspon- 
dence, we shall demonstrate in this communication that 
the exact quantum Neumann eigenstates can be used 
to expand and evolve an arbitrary classical projected- 
density. 

The paper is organized as follows. In section|TTl we deal 
with the projected density evolution operator for billiards 
and discuss some of its general properties. In section UIIl 
we explicitly show how a projected-density evolves in an 
integrable system and summarize the results Q on the 
eigenvalues and eigenfunctions of £p for billiards in gen- 
eral. Section IIVI contains numerical demonstration that 
quantum states can be used to evolve classical projected 
densities in an integrable as well as a chaotic system. Fi- 
nally a summary and discussion of the results is provided 
in the concluding section. 



boundary. Depending on the shape of the boundary, bil- 
liards exhibit the entire range of behaviour observed in 
other dynamical systems. The quantum billiard problem 
consists of determining the eigenvalues and eigenfunc- 
tions of the Helmholtz equation 



V^ij{q) + k^ijiq) = 



(11) 



with 4'{q) — (Dirichlet) or n.Wtp ~ (Neumann 
boundary condition; n is the unit normal) on the bound- 
ary. Apart from being a paradigm in the field of classical 
and quantum chaos, billiards have relevance in a vari- 
ety of contexts including the motion of ultra-cold atoms 
confined by laser beams (the so called "optical billiards" ) 
[13, [Hi- The Helmholtz equation describing the quantum 
billiard problem also describes acoustic waves, modes in 
microwave cavities and has relevance in studies on "quan- 
tum wells" , "quantum corrals" , mesoscopic systems and 
nanostructured materials. The Neumann boundary con- 
dition has important manifestation in acoustic waves, 
surface water waves, TE modes in cavities and modes 
of a drum with stress-free boundaries [l^] as well as ex- 
citations of Bose- Einstein condensates in billiards • 

As the motion inside the enclosure is free in a billiard, 
the magnitude of the momentum is conserved. The kernel 
for a billiard is thus 

Kp{q,q',t) = ^ J dpdp'6{q~q'')5{p~p'')S{pl-p'') 

] * t ^^^^ 

where q' ~ q' {q',p'), p' — p' {q',p') is the position 

after time t of an initial phase space point {q',p')- It is 

assumed henceforth that the kernel Kp depends on the 

parameter Eq — Pq. 

It is convenient to treat p in polar coordinates {p, ip) 

where p is the magnitude and ip is the angle that p 

makes with a given axis. Transforming from {px,Py) to 

{p, ip) and noting that p is conserved, the p' integration 

in Eq. can be simplified as (see appendix E| : 



dp'Jp'y5{p.-p'J{q',p'))5{py-p'y\q',p'))h{q',p') 
dip' S{p~p'')h{q',^';p) = h{q',^*;p) (13) 



where p'J^{q',p') and p'y(q',p') are the x and y compo- 
nents of the momentum at time t for the initial phase 
space coordinate (q' ,p') and p* is the value of p' for 
which p'*- {q' , p' , p) = pi [1J|. In the above. 



II. THE PROJECTED DENSITY EVOLUTION 
OPERATOR FOR BILLIARDS 

We shall henceforth restrict ourselves to billiards where 
the eigenfunctions of C^p have interesting properties. In 
a classical billiard, a point particle moves freely on a 
plane inside an enclosure and reflects specularly from the 



h{q', p-p) = 5{q - q'*{q', p,p))S{p' - pi). (14) 
The kernel thus simplifies as 

Kp{q,q',t) = ^J pdpdp, 5{q-q'\q'-p,p*))5{p^ -pi). 

(15) 
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On carrying out the p integration, the kernel can finally 
be expressed as 



Kp{q,q',t;po) 



2tt 



dip 5{q~q'\q';pQ,ip*{p))). 



(16) 

Thus, the time evolution of the projected-density can be 
determined using Eq. ^ and Eq. (fT6|) . 

We list below some properties of the projected-density 
evolution operator, £p: 

1. £p preserves positivity of the density. To see this, 
note that the delta function in the kernel {Kp) 
picks up values from the initial density, which, is 
positive everywhere by definition. Thus the inte- 
grand remains positive at all times. 

2. The operator £p is not multiplicative i.e. /Ip"*"*^ ^ 
^p" o r*^ . A simple way of understanding this 
is by noting that evolution of the projected den- 
sity is such that trajectories move out isotropi- 
cally from each point of the initial density. Thus, 
£p o (£p o p{q)) evolves the initial density isotropi- 
cally to £p o p[q) followed by an isotropic evolution 
of this density for a time i2- This is different from 
an isotropic evolution of the initial density for time 
ti -\- 12 (see appendix [B|) . 

A consequence of nonmultiplicative evolution is 
that the eigenvalues of £p are not of the form e'*''*. 

3. The eigenvalues and eigenfunctions of £p have in- 
teresting properties in billiards [1, Q as elaborated 
in the following section. It is found that the eigen- 
functions are (at least) well approximated by the 
quantum Neumann eigenfunctions while the eigen- 
values are {/(\/i?„wi)} where {En} are well ap- 
proximated by the quantum Neumann eigenvalues. 
Here v is the speed of the classical particle and f[x) 
is the asymptotic form of the Bessel function Jo (a;). 



III. EIGENVALUES AND EIGENFUNCTIONS 

OF C'p 



CM') - / dOo dIop{Io,eo)Pr^Mo,0o) 

d0o p(I',0o)e-™-^«. (18) 



//2 



(2^) 



The evolution of p{I, 9) in terms of the eigenvalues and 
eigenfunctions of £* can thus be expressed as 



C'op{I,e)^Yl j dI'Cr,{I') Pi'ALO) e™"(^')* 

n •' 

(19) 

Using Eq. ([3]) for the eigenfunctions, one can express 
Eq. HID as 



C'op{I,6) = J2 I dI'Cr,{I')5{I-I')~pr,e'^-^^''^' 

n 

= ^C„(/)e'"'^(-f)*p„(0) (20) 



where 



dl pi'^rv = - — "77^. 

(27r)^/^ 



(21) 



Our interest here is in initial densities that are con- 
strained to the constant energy surface having the form 



PE{I,e)^S{E-H{I))g{e). 
Thus using Eqns. HH]), ^ and 



(22) 



C*opEiI,e) = E(/ d9o g{eo)p*jeo) 

X 6{E~ H{I))e'"-'^^'^*pr,{e) (23) 

Finally, in order to project pE to 6 space, we must inte- 
grate £* o pe{I, d) over I to obtain 



A. Integrable Case 



For integrable systems, the eigenmodes of £p can be 
determined easily. In terms of action-angle variables 
{/, 0}, an arbitrary initial density p{I,6) can be ex- 
panded as 



p(J,0)=^ J dJ'C„(/')P/^r^(/,0) 



(17) 



where pi\n are the eigenfunctions of £* (Eq. [3]) and the 
coefficients 



CpOppiO) - Ppi0o)p*M'^ Pr^iO) 



dIS{E - H{I))e 



in.uj{I)t 



^(7„ A^ Pr,{9) 



(24) 
(25) 



where pp{6) = J dl pe{I, S) = pg{d) and the quantities 
within (...) and [. . . ] are respectively the coefhcient of 
expansion ((?„) of the projected density pp{0) in terms 
of Pn and the corresponding eigenvalue (A^) of the pro- 
jected evolution operator £p. 
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To illustrate this, we shall consider the rectangular bil- 
liard which is an integrable system. The hamiltonian 
expressed in terms of the actions, h-h is -ff(/i,/2) = 
7r^(/^/L^-|-/|/L2) where ii, L2 are the lengths of the two 
sides and the mass ra = 1/2. With Ii = \/ELi cos{ip)/iT 
and I2 = ^pEL^ sin((^) /tt, at a given energy, i?, each torus 
is parameterised by a particular value of (p. Transforming 
from Ji, I1 to i?, we have 



PP 



m 

L1L2 



L1L2 



27r2 



J dE dip d{Eo ~ E) 



g{e) = figie). 



(26) 



The eigenvalues of £p (denoted as A^) are 



L1L2 



dEdipS{Eo - E)e 



27Tit-\/E[^^ cos sin ip] 



= JoiV Enl) 



(27) 



where En = 'n'^{n\/L\ + "-2/^2)' tan/i„ = n^lrii, I = 
and Jo is the zeroth order Bessel function. Thus, from 
Eq. [^H it follows on putting ppifi) = Pm{d) that 



C^P o PmiO) ^ Joi\/E^l) Pm{d) 



(28) 



as Cn — 5m.n- Note that m = (toi,TO2) are natural 
numbers so that e^^"^^^^^™^^^^ are also the quantum Neu- 
mann eigenfunctions in 9 space while Em are the quan- 
tum Neumann energy eigenvalues. 

We end this sub-section with an observation. For the 
integrable case, the density 



Pe,^ 



1 in.e 

-^(^o-i^(I))(^ 



(29) 



on the constant energy surface is an eigenfunction of £* 
as £* o pe^ti — e'^-^V-B.Ti- Since 



J dl S{Eo-H{I))= 1^, 
the momentum projection of pE,n is 



(30) 



An.Q 



dl PE,: 



Pr.{e). 



(27r)//2 

Note that Pn{d) is an eigenfunction of £p i.e. 



(31) 



C'p O p„(0) = Jdl Co pE,r. = A^ p„(0) (32) 



In contrast, for an eigenfunction on the torus (Eq. [3]), 
momentum projection yields 

j dl to pi,^^ = p^{e) j dI6{I - J')e'"'^(^)* 

= e™-'*'(-^')*p„(0) A^e*"(-^')* (33) 

Thus the properties of the projected operator depend on 
the function space to which the phase space density be- 
longs. 

B. Billiards - The Quantum Connection 

It is instructive to see that the quantum Neumann 
eigenfunctions in q space are also eigenfunctions of £p. 
As an illustrative example, consider a particle in a one- 
dimensional box with walls placed at g = and q = L. 
To find the evolution of the quantum Neumann eigen- 
fmiction V'ra(?) = e.'^^"'' + e"''^"'^, kn = mr/L under £p, 
note that in this 1-dimensional case there are only two 
possible values of (p i.e. = 0, tt. Thus 



[£*(+) + £*(-)] 



(34) 



where £*(±) refer to positive [p — 0) and negative {p — 
it) velocity. The classical evolution for positive velocity, 
o ipn{q), is given by 



(35) 



where q *{+v) is the position at time —t with initial 
position q and initial velocity +v. Similarly, C*{—)oipn{q) 



(36) 



with the — sign in £*(— ) denoting negative velocity. Note 
that the flow is such that the velocity changes sign at 
every reflection from the walls placed at g = and q = L. 
For the flow q~*{+v), the reflections occur at t+ = (q + 
nL)/v so that for < t < q~*{+v) = v{t - t^) = 
vt — q. Similarly, for the flow q^'^{—v), the reflections 
occur at ~ {L — q + nL)/v and for < t < t\ 
q-\~v) ^ L ~ v{t - 
that 



2L — vt — q. It follows hence 



£*(±) o V'„(g) = e*'="(-« ± + e-''^"^-* ± (37) 
for all t. Thus 
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£*P o ipn{q) = C0ii{knVt)ipn{q). 



(38) 



In other words, the quantum eigenfunction is also an 
eigenfunction of the projected classical evolution oper- 
ator £p. 

For general 2-diniensional billiards, there are two ap- 
proaches for determining the eigenmodes of £p [1, 
Both rely on polygonalization of the billiard boundary 
[isl ITg} and the argument that as the number of sides of 
the polygon is increased, the modes of the smooth bil- 
liard are approximated better. In the first approach (17j , 
the trace of the £p is related to the trace of the energy 
dependent quantum Neumann propagator while in the 
second, a plane wave expansion enables one to conclude 
that there exists a correspondence between the quan- 
tum Neumann eigenmodes and the modes of £p [11] . 
We merely state here the result and refer the interested 
reader to [1 Q. For i > 0, 



>Cp o V'n(q) = f{k„vt) i)n{q) 



(39) 



where En — are well approximated by the quantum 
Neumann energy eigenvalues, f{x) = -^/2/7nz: cos(x— 7r/4) 
is the asymptotic form of Jo{x) and ipnil) are well ap- 
proximated by the quantum Neumann eigenfunctions. 
For some integrable systems, it is possible to show that 
{En} and {"(An} are the exact quantum Neumann eigen- 
values and eigenfunctions. 



IV. EVOLVING CLASSICAL DENSITIES WITH 
QUANTUM STATES: NUMERICAL EVIDENCE 

As stated earlier, the purpose of this communication 
is to demonstrate that quantum Neumann energy eigen- 
functions {tpniq)} can be used to evolve an arbitrary 
classical momentum- projected density pp{q), despite the 
fact that Eq. is generally approximate. 

The completeness of the quantum Neumann eigenfunc- 
tions allows us to expand pp{q) as: 



Cn'ipniq) 



(40) 



where 



C„ = y dq pp{q)i)l,{q). 



(41) 



In the above, it is assumed that / dq ?A^(g)V-'Tn('Z) = 
5m,n- On using Eq. ([M)) . we have 



£*P o pp{q) = CMVE^vt)ij;^{q). (42) 



In order to account for the behaviour at small t, we use 
the Bessel function Jq since f{x) is an asymptotic form of 
Jo (a;). Note that the evolution of the projected-density, 
pp{q), under Eq. ([^^ should reflect isotropy in momen- 
tum. 

To demonstrate that Eq. indeed gives us the cor- 
rect classical evolution, we choose points (initial condi- 
tions) in <7-space distributed according to the initial den- 
sity pp{q) and evolve them in time using classical equa- 
tions of motion. The evolution of these points give us the 
density pp{q, t) at any time t. This is compared with the 
density obtained from Eq. ([42]) . 

In both the examples considered below, we shall con- 
sider a hat function as the initial density for convenience. 
The function takes the value unity within a square strip 
of side AL centred at a point q and is zero outside. Thus, 
all initial conditions lie uniformly within the square strip. 
For convenience, we choose the velocity v = 1. 

As a first example, consider a particle in a rectangular 
box of side ii = L2 = 1-0. Fig. ([T^) shows the initial 
density centred at (0.2,0.2) with AL = 0.15. The first 
5000 quantum Neumann states have been used to expand 
the initial density. Figs. (Pd) to ([TJ3) are the densities 
at time t = 0.25,0.5,0.75,1.0 respectively as computed 
using Eq. (|42)l . For comparison, Fig. [2] is a 2-dimensional 
view of the density evolved using classical trajectories 
with their momentum distributed isotropically Pi. At 
time t = 0, the points lie in a square patch as shown 
in Fig. ([2^) while Figs, ^p) to ([2};) show the points as 
they evolve from Fig. (HJi) at times t = 0.25, 0.5, 0.75, 1.0 
respectively. Note that the evolution in Fig. 1 reflects 
isotropy in momentum as expected. 

We next consider the chaotic Stadium billiard consist- 
ing of two parallel straight segments of length 2 joined 
on either end by a semicircle of unit radius. The ini- 
tial density is again centred at (0.2,0.2) with AL = 0.15. 
In this case, the quantum Neumann eigenfunctions have 
been computed numerically using the boundary integral 
method. We use the first 1000 eigenfunctions to expand 
the initial density and evolve it using Eq. ([i^ . Fig. [3] 
is similar to Fig. [l]but for times t=0,l,2,3,4 and 5 while 
Fig. [4] (similar to Fig. [2]) shows the classical evolution of 
the initial points at these times. The similarity between 
the two evolutions is again obvious despite the approxi- 
mate nature of the eigenfunctions and eigenvalues used. 

At longer times, both evolutions lead towards a uni- 
form distribution in q space. Fig [5] shows the density 
evolved using Eq. [55] at times t=10, 20 and 30 corre- 
sponding on an average to 4.5, 9.0 and 13.6 bounces. 
This is compared to the classical evolution of the initial 
points as these times in Fig. [6l Note that the final in- 
variant density should assume a value p ~ 0.00315 for 
the initial density considered. Clearly, the density evo- 
lution using Eq. 221 approaches this value with fluctua- 
tions around it that diminish with time (see Fig. [5])- At 
i = 30, the density appears to be nearly uniform in both 
cases (Fig. ([St) and Fig. ([5J:)). The decay to the cor- 
rect (uniform) invariant density is not surprising as the 
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Neumann ground state eigenfunction is uniform and the 
corresponding eigenvalue i5o = so that Jo{y/Eovt) — 1. 
However, the finer structures are harder to discern at 
fonger times. 

A closer inspection of the classical evolution however 
shows that that the density is nonuniform even at t = 
30. Fig [3 shows the evolution at t = 20, 40, 60 and 110 
using a finer representation (dot size) of a point. In the 
first two cases, the departure from uniformity seems ev- 
ident while at t = 60 (average of 27 bounces), a closer 
inspection shows the presence of four patches of higher 
density. Finally at t = 110, the density appears to be 
uniform corresponding to an average of 50 bounces per 
trajectory. The classical evolution of a Gaussian pro- 
jected density centred a,t q — (0.2, 0.2) is shown in Fig. [H] 
at t = 0, 20, 40 and 60. Note that at t = 40, the density 
is non- uniform while at t = 60 (average of 27 bounces), 
the density appears to be uniform. The decay rate thus 
depends on the initial projected density. 

There are two sources of errors in the stadium bil- 
liard evolution. In the first place, the quantum Neumann 
eigenfunctions are approximate classical eigenfunctions, 
and, the quantum Neumann eigenvalues {En} used in 
fiVE^tv) are close to but not exactly equal to the clas- 
sical values. The latter is borne out by the fact that the 
position of peaks in the fourier transform of Kp are close 
to but not exactly at the quantum Neumann eigenvalues 
[13] ■ The second source of error arises from the the trun- 
cation of the basis. As associated problem is the fact that 
the quantum Neumann eigenfunctions used in the study 
have been evaluated numerically and hence have errors. 
In the stadium billiard a thousand eigenstates have been 
used to expand the initial density while in the rectan- 
gle (where an analytic expression is used) five thousand 
eigenstates have been used to expand and evolve the den- 
sity. The oscillatory background of the initial density in 
Fig is a consequence of truncation error and the er- 
ror in numerical evaluation of the exact Neumann eigen- 
functions. This, coupled with the first source of error 
mentioned above, can even lead to negative values of the 
density at some points. 

Despite these problems, evolution using the quan- 
tum Neumann eigenfunctions and the eigenvalues 
{f{\/Entv)} captures the classical trajectory evolution 
fairly accurately. The finer structures are reproduced 
well for about 5 bounces while at longer times, evolu- 
tion using the quantum Neumann eigenstates leads to 
the correct invariant density. 



V. DISCUSSION AND CONCLUSIONS 

We have dealt with the evolution within a billiard en- 
closure of an initial classical density on the energy shell 
that is isotropic in momentum. We have constructed 
an appropriate classical evolution operator (£p) for the 
momentum projected density and shown that its approx- 
imate (sometimes exact) eigenvalues and eigenfunctions 



have a one-to-one correspondence with the quantum Neu- 
mann eigenvalues and eigenfunctions. Based on this, we 
have demonstrated that for both the rectangular and sta- 
dium billiards, expansion of the initial projected-density 
using the quantum Neumann eigenfunctions as eigen- 
functions of £p, captures the classical evolution fairly 
accurately. While the finer structures are reproduced at 
shorter times (about 5 bounces for the stadium billiard) , 
evolution using the quantum Neumann eigenstates leads 
to the correct invariant density at longer times. Note 
that the spectrum of the evolution operator £p does not 
fix a time scale in which the initial density decays to the 
invariant density. 

It is instructive to discuss the difference between evolu- 
tion due to C*p and a purely quantum evolution. Specifi- 
cally, we may ask whether the hat function that we con- 
sidered as the initial classical density, can be evolved 
quantum mechanically in an identical fashion at least for 
short times. The initial quantum state ip{q, t — 0) in this 
case is again a hat function (as it is the square root of 
the initial density, p{q) ) and its evolution is governed by 

ijiq,t) = Y,Cne'''"'MQ.O) (43) 

n 

where c„ — Vri(9, 0)c?q. Here A refers to the region 
in configuration space where the initial wavefunction is 
non-zero. Eg. 1431 however contains no information about 
the energy or momentum of the underlying classical tra- 
jectories. Thus, evolution due to Eq. 23] cannot approx- 
imate classical evolution. At the semiclassical level, this 
is due to the fact that trajectories with all possible en- 
ergy contribute to the propagator and a correspondence 
with classical evolution is possible only when the Wigner 
transform of the initial wavepacket is localized in phase 
space. In contrast, evolution due to >Cp has information 
about the magnitude of the momentum {v in the eigen- 
value) and the isotropy in momentum is built into the 
kernel Kp and reflects in the form of the eigenvalue (the 
Bessel function Jq). 

Another important distinction concerns the decay 
rates of initial phase space densities that are direction- 
ally (or if) dependent and independent (isotropic or 
uniform in if). To explore this further, we shall use 
the operator V to denote momentum projection. For 
the anisotropic or directionally dependent initial density 

{VoC')p{q,p) = VY,CnC'Mq,p) 

n 

= ^C„e""*7'0„(q,p) 

n 

= 5]C„e""*^„(q). (44) 

n 

Here 0„(q) = V(j)n{q,p), <j)n{q,p) are the eigenfunctions 
of £p and the coefficients C„ can be determined using the 
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eigenf unctions of C\ the adjoint of C If the leading a„ 
has a negative real part as in case of hyperbolic systems, 
the decay to the invariant density is exponential. 

In contrast, an isotropic initial density such as 
p{QtP) = p{QtP) — 9{q)^{Eo ^ P^) cannot be expanded 
in eigenfunctions of £*. To understand this, note that 
the only eigenfunction of £* which is uniform in ip (i.e. 
isotropic) is the invariant density which is uniform in 
both q and (p. All other densities that are initially uni- 
form in (f, lose their uniformity (in ip) as they evolve with 
£*. In other words, £* does not possess any eigenfunc- 
tion that is uniform in tp but not uniform in q. Thus 
g{q)S{Eo — p^) cannot be expanded in an eigenbasis of 
£*. The decay to the invariant density must therefore be 
dictated by the eigenvalues of £p since 



{VC) g{q)5{E^^p^)^C'p pp{q) 



(45) 



where pp{q) — J dp p{q,p) — T^giq)- Thus the decay 
rates of initial densities that are isotropic must differ from 
the decay rates of anisotropic densities. 

A few conclusions based on the results of the preceding 
sections are listed below: 

• The evolution of a momentum projected density 
using the (approximate) eigenstates of the momen- 
tum projected evolution operator Cp captures the 
classical evolution fairly accurately. The approx- 
imate eigenfunctions are the quantum Neumann 
eigenfunctions while the eigenvalues are related to 
the quantum Neumann eigenvalues. 



• The evolution operator Cp is non-multiplicative i.e. 
r^^+ti ^ Qt2 o £ti^ Tj^^g^ ^j^g eigenvalues do not 

have the form e'*'**. 



• £p preserves positivity of the initial density. How- 
ever numerical errors and the approximate nature 
of the eigenvalues and eigenfunctions used in the 
evolution (Eq. W2^ can lead to negative values at 
some points. 

• As the small t behaviour of the density is repro- 
duced by Eq. (|42|) . the correct form of the approx- 
imate eigenvalues of £p is { Jo(\/ Envt)\ where En 
are the quantum Neumann energy eigenvalues. For 
integrable billiards such as the rectangle, these are 
the exact eigenvalues. 

• For non-polygonal billiards such as the stadium, 
proof of the existence of a correspondence was in- 
direct, based on the limit of a large number of sides 
(of the corresponding polygonalized billiard), and, 
numerical evidence based on a few eigenstates 

The results presented here show that this corre- 
spondence must hold for all states as the evolution 
of the density using the quantum Neumann eigen- 
states faithfully follows the trajectory picture. This 
puts the correspondence on a firmer footing. 



The relaxation of an arbitrary projected-density to 
the uniform, steady-state density can be predicted 
reasonably using the approximate eigenvalues and 
eigenfunctions of the projected Perron-Frobenius 
operator, £p. 



APPENDIX A 



We shall show here that in polar coordinates (p, ip) , 



1 S{p-p'')h{q',p')dp' = J V d{ip~ip"Mq',p,p'). 

(Al) 

where = ip*{q' ,p, tp') is the polar angle of the momen- 
tum vector p' at time t for an initial phase space point 
{q' ,p, p>'). Suppressing q', the integral can be expressed 
as 



S{p cos ip — p' cos ip'*)S{psm (p ~ p' sin tp'^)h{p')dp' 

(A2) 

where it is assumed that for billiards p'* — p' . Obviously 
the only value oip' that contributes to the integral is p' ~ 
p and (/s'* = Lp. In order to carry out the p' integration, 
we may rewrite (IA2[) as 



5(p^^ -p') 
^"'^'^ „ ' S{psmp^-p'sm^")h{p')dp' (A3) 

cos <y9" 



On performing the p' integration, we have 



7 / • lf\ pcosip , , 
dip d [p sm ip — p cos ip tamp ) ^ — -n(p,ip ) 



cos^ p>" 



, 5{t&nip- tan ip'*-) 
dip 5— h{p,ip) 



COS^ ip'* 



dip' Sip-p=")hip,p,'). 



(A4) 



APPENDIX B: NON-MULTIPLICATIVE 
PROPERTY OF Cp 

The evolution of momentum projected isotropic initial 
densities on the energy shell is governed by £p : 



Cp o pp{q) ^ J dq' Kp{q, q' , t) pp{q'). (Bl) 
The non-multiplicative property of £p implies 

j dq" Kpiq, q", h) Kp{q" , q', h) ^ Kp{q, q' , h + h) 

(B2) 
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This is easiest to visualize when the points q', q and 
the times ti and ^2 are such that such that there is no 
encounter (reflection) with the boundary of the bilhard. 
As the momentum vector does not change direction, (p* — 
(fi in Eq. [161 Using Eq. [T6]for the projected kernel Kp, 

J dq" Kp{q,q",t2) Kp{q" ,q' M) 

= / d^d^'dq" S{q-f,'^q",po,v')) 

^S{q" - f,'^{q',P0,^)) (B3) 
= ^ / d^d^' Siq~ f,'^f,'-{q',ip),^'))iB4) 

^ l-J d^5(q-/*i+*^(q',po,^)) (B5) 

where fg (also denoted by q* in the text) is the q- 
component of the flow /* . 

The last step can be best understood in terms of Fig. [9] 



which illustrates the difference in the two time evolutions. 
The first is a two step evolution from q' to q" in time ti 
followed by an evolution from q" to q in time t2- These 
are indicated by solid lines terminated by arrows. The 
kernel determines the angles ip and tp' which connect q' 
with q via the intermediate point q" . Depending on the 
angle tp, the point q" can lie anywhere on the inner circle 
(Ci). Any point q on the circle C2 (with centre on Ci) 
has a contribution at time ti + 12 from the initial density 
at q' . The second evolution is for a time ti + t2 starting 
from q' . The kernel for this evolution contributes when 
q lies on the outer circle (C3). 

Thus, £p o £p ^ £p"'"*^. A consequence of the non- 
multiplicative nature of >Cp is that its eigenvalues cannot 
be of the form e'^'* i.e. the eigenvalues of £* are dis- 
tinct from the eigenvalues of £p when the initial phase 
space density is isotropic and constrained to the constant 
energy surface. 



[1] P. Cvitanovic et al. Chaos: classical and quantum, 

Jittp:/ /chaosbook.org 
[2] A. Lasota and M. MacKey, Chaos, Fractals, and Noise; 

Stochastic Aspects of Dynamics, Springer- Verlag, Berlin, 

1994. 

[3] P. Gaspard, Chaos, scattering and statistical Mechanics, 

Cambridge University Press, 1998. 
[4] J. Wilkie and P. Brumer, Phys. Rev. A 55, 29 (1997). 
[5] D. Biswas, in Nonlinear Dynamics and Computational 

Physics, ed. V. B. Sheorey, Narosa, New Delhi, 1999; 

chao-dyn/9804013 
[6] W. T. Lu, W. Zeng and S. Sridhar, Phys. Rev. E 73, 

046201 (2006). 

[7] These are not the same as the eigenvalues and eigen- 
functions of See section HTl and appendix iBl for non- 
multiplicative property of £p, section Hill for an illustra- 
tion using integrable billiards and the discussion in the 
concluding section. 
[8] D. Biswas, Phys. Rev. Lett. 93, 204102 (2004). 
[9] D. Biswas, Pramana - J. Phys, 64, 563 (2005). 

[10] V. Milner, J.L. Hanssen, W. C. CampbeU and 
M. G. Raizen,Phys. Rev. Lett. 86, 1514 (2001). 

[11] N. Friedman, A. Kaplan, D. Carasso, N. Davidson, Phys. 
Rev. Lett. 86, 1518 (2001). 

[12] T. A. DriscoU and H. P. W. Gottlieb, Phys. Rev. E 68, 
016702 (2003); H.-J. Stockmann, Quantum Chaos - an 
Introduction, Cambridge University Press, Cambridge, 
1999. 

[13] C. Zhang, J. Liu, M. G. Raizen and Q. Niu, Phys. Rev. 
Lett. 93, 074101 (2004). 

[14] In polygonalized billiards, if'^{q' ,Lp' ,p) — Lp' for a tra- 
jectory in the unfolded space which is produced by suc- 
cessive reflections of the billiard domain about the edge 
that the trajectory strikes. Thus in the unfolded space, a 
trajectory is a straight line. 

[15] J. L. Vega, T. Uzer and J. Ford, Phys. Rev. E 48, 3414 
(1993). 

[16] D. Biswas, Phys. Rev. E 61, 5073 (2000). 



[17] D. Biswas, Phys. Rev. E 63, 016213 (2001). 

[18] D. Biswas, Phys. Rev. E 67, 026208 (2003). 

[19] For a better comparision, for each initial point in q, an 
isotropic distribution of momentum should be consid- 
ered. 



9 





FIG. 1: The evolution of a projeeted-density, p{q), in the rectangular billiard is computed using the exact quantum Neumann 
eigcnfunctions as the cigcnfunctions of the projected Perron- Frobenius operator, £\,. The initial density is shown in (a) while 
(b) to (e) show the density at t=0.5, 1.0, 1.5, 2.0 respectively. The evolution reflects isotropy in momentum. 




FIG. 2: The initial points chosen according to the density in Fig. la is shown in Fig. 2a. These are evolved classically using 
momentum direction (if) distributed uniformly in (0, 2n). The distribution of points at time (b) 0.5 (c) 1.0 (d) 1.5 and (e) 2.0 
can be compared with the densities in Fig. lb to le respectively. 




FIG. 4: As in Fig. 2 for the stadium billiard at times (a) t=0 (b) t=l (c^ 



■) t=2 (d) t=3 (e) t=4 and (f) t=5. 
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FIG. 6: The initial density of Fig. 4a evolved classically using momentum directions distributed uniformly at times (a) t=10 
(b) t=20 and (c) t=30. 




FIG. 7: The initial density of Fig. 4a evolved classically using momentum directions distributed uniformly at times t=20 
(top-left), t=40 (top-right), t=60 (bottom-left) and (bottom-right) t=110(bottom-right) . 




FIG. 9: Comparison of two time evolutions in a billiard enclosure (boundary not shown). 



